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Steady but generally unstable solutions of the 2D Boussinesq equations are obtained for no-slip 
boundary conditions and Prandtl number 7. The primary solution that bifurcates from the con¬ 
duction state at Rayleigh number Ra ~ 1708 has been calculated up to Ra « 5.10® and its Nusselt 
number is Nu ~ 0.143 with a delicate spiral structure in the temperature field. Another 

solution that maximizes Nu over the horizontal wavenumber has been calculated up to Ra = 10® 
and scales as Nu ~ 0.115 for 10^ < Ra < 10®, quite similar to 3D turbulent data that show 

Nu - 0.105 i?a0'3i in that range. The optimum solution is a simple yet multi-scale coherent solu¬ 
tion whose horizontal wavenumber scales as 0.133 That solution is unstable to larger scale 

perturbations and in particular to mean shear flows, yet it appears to be relevant as a backbone for 
turbulent solutions, possibly setting the scale, strength and spacing of elemental plumes. 
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Rayleigh-Benard convection is the buoyancy-driven motion of a fluid contained between two horizontal 
plates and heated from below. It is a paradigmatic problem with a rich array of nonlinear physics from 
deterministic chao^ and period doubling to pattern formatiorP and turbulenc^. A fundamental issue is 
to determine how the heat flux R scales with the temperature difference AT between the bottom and top 
plates. There is a critical NT,, such that heat transport is by conduction with R ^ AT for AT < AT^ but 
macroscopic fluid motion develops for AT > ATc and enhances heat transport. Bifurcations take place as 
AT is increased, leading to turbulent flow and R ~ (AT)^/® and perhaps even as high as R ^ (AT)^/^. 
Understanding these bifurcations and the actual asymptotic scaling of R have been the subjects of many 
experimental, theoretical and numerical studies focusing on turbulent convection. Here, we consider coherent 
convection - steady flows that may be the permanent form of the coherent structures (‘plumes ’, ‘thermals ’ 
and large scale ‘winds ’) permeating turbulent convection. 

In Rayleigh-Benard convection, the fluid is contained between two infinite horizontal plates and the 
Boussinesq approximation is made so the governing equations are 

dv 

— + {v ■'V)v+ 'Vp = ga^Ty + (1) 

BT 

— + (u • V)T = «V®T, (2) 

for the velocity v and the temperature T. Incompressibility V • u = 0 is maintained by the kinematic 
pressure p and —gy is the constant acceleration of gravity. The fluid density is p « po (1 ~ with 

c«v — 0 the volumetric thermal expansion coefficient and po a constant reference density, ^' > 0 is the 
kinematic viscosity and k > 0 the thermal diffusivity. The bottom plate temperature is fixed at AT/2 and 
the upper plate is at —AT/2. The conduction state is the solution d = 0, T = —yAT/H to Q and ([^, 
with —h < y < h = H/2. This is a steady solution for all values of the parameters, but Rayleigh showed 
that it is linearly unstable when the Rayleigh number Ra = ga^,AT H^/{vk) is larger than a critical value 
Rac, independent of the Prandtl number Pr = v j k. 

We consider no-slip boundary conditions for which Roc ~ 1708 for horizontal wavenumbei!^ Oc = 3.116/77. 
Convection develops and enhances heat transport for Ra > Rac- For statistically steady solutions, integra¬ 
tion of ([^ yields the heat flux,l^ 

= RqP{vT) (3) 

V=±h 

where Rq = kAT/H is the heat flux in the absence of fluid motion and v = y -v is the vertical velocity with 
T and (*) denoting horizontal and domain averages, respectively. 


R = -k 


dT 

dy 
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A classic scaling argument is that heat transport is determined by marginal stability of the thermal bound¬ 
ary layers and of the mean temperature gradient in the interiorPA simple version of this argument assumes 
an isothermal T = 0 interior with boundary layers of thickness S such thalP = ga^(AT/2)6^/ii'K) « 
1708/16, yielding heat flux 

n ~ \ ' (At)4/3 0.084 (4) 

6 \ 1708 V) d 

independent of Pr, where the Nusselt number Nu = 'H/'Hq. This scaling also follows from assuming that 
the heat flux becomes independent of the height H as Ra —>■ oo. 

Another simple argument imagines a fluid at temperature T = 0, with cold plumes at temperature — AT/2 
free-falling a distance H = 2h from top to bottom plates (and hot AT/2 plumes ‘free-rising’ from bottom 
to top) at average speed V = y/g'h, where g' = ga^AT/2 is the reduced gravity. This yields the heat flux 

H ~ VAT/2 = (AT)3/2 Nu ^ ^ = ^{Ra . (5) 

This inertial scaling also follows from the standard turbulence assumptiorP that heat transport becomes 
independent of v and k in the limit Vh/v —> oo, Vhj k —> oo. 

Kraichnan’s mixing length theorjP yields various scalings in distinct regions of the (i?a, Pr) parameter 
space. He obtains Nu ~ Ra^^^, for fixed Pr, but predicts a transition to Nu ~ Ra}^^{\n Ra)~^/'^ for very 
large Ra > 10^^ when the shear boundary layers would be turbulent. Grossmann and Lohse’s comprehensive 
scaling theorjffSI incorporates classic Kolmogorov scaling of velocity and temperature dissipation rates in the 
bulk together with thermal and viscous boundary layers. Their theory fits many experimental data sets.^ 
They predict a transition to an ‘ultimate’ regime similar to Kraichnan’s for Ra > 10^^ and argue that their 
log correction would yield an effective scalinj^^ of about Nu ^ Ra^-^^ in the range 1 0^^ < Ra < 10^®. 
Some experimental results show transition to an ultimate regim^i^lti^ while others do not. b'^l^s l However, all 
the data is well-fitted by Nu ~ 0.105 for Ra < 10^^ and that is the regime for which we report on 

unstable coherent states, in particular optimum transport solutions with Nu ~ 0.115 

Rigorous upper bounds on heat transport for no-slip boundary condition^iSH^S] yield Nu — 1 < 0.026 Ra^/"^ 
as Ra —>■ o o, for any Pr, showing that the free-fall scaling ([^ cannot hold for large Pr. Rigorous bounds 
for free-slijP^HU yield Nu — 1 < 0.106 that conflicts with (§ for any Pr. 

We consider 2D flow with v = u{x,y,t)x + v{x,y,t)y and use h = H/2, V = ^/g'h, r = h/V, AT/2 as 
our characteristic length, velocity, time and temperature scales, respectively. Eliminating p by taking the y 
component of the curl of the curl of Q yields 


dtV^v = V + dlT + dx {vN^u - uV\), 

(6) 

dtT = K V^T - (u • V)T, 

(7) 

dtu = V dyU — dyWv 

(8) 


where v and k are now non-dimensionalized hy Vh = ^g'lp and relate to the Rayleigh and Prandtl numbers 
as 


r = 



4 

yjRa Pr 


(9) 


All results in this paper are for Pr = 7 (water). The horizontal velocity u is obtained from v and dxU+dyV = 0 
except for its horizontal average u{y,t) that is determined by ([^, where the overbar denotes a horizontal 
average. Equations (§, 0, 0 are considered with no-slip boundary conditions 

u = 0, dyV = 0, T = =f 1 at y = ±1, (10) 


together with periodicity of period T = 27r/a in the x direction. 

We look for steady solutions, dt = 0, that bifurcate from the conduction state at Ra « 1708 for wavenum¬ 
ber a « 3.116/2 corresponding to aspect ratio L/Tl « 2.016. Those convective solutions obey mirror 
symmetry 


[u,v,T]{x,y) = [-u,v,T]{-x,y) 


( 11 ) 
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FIG. 1. Nu vs. Ra for steady primary branch (PB) with L/H = 2 and optimum branch (OB) bifurcating from 
conduction state at Ra = 1708, Nu = 1. Primary branch bifurcates to a time-periodic solution near Ra = 53 000, 
max and min Nu achieved by the periodic solution are plotted. Unstable steady state is dashed. 


as well as shift-reflect symmetry 


[m, n, T]{x, y) = [u, -v, -T]{x + ^,-y)- 


( 12 ) 


Equations @,0 are discretized using a Fourier expansion in x and Chebyshev integratiorP^ in y with time- 
marching to steady states. Two distinct codes have been written, code 1 uses a 3rd order time-accurate 
scheme,^ code 2 uses a semi-implicit forward-backward Euler time discretization together with Chebyshev 
tau of the 2nd kincP^ for the v equation. Code 2 typically uses inconsistent time integration, with smaller 
time steps for smaller wavenumbers and for the v equation than for the T equation, to speed up or enable 
convergence to steady state. Both codes are dealiased in x and y using the 2/3 rule. The results of both 
codes overlap or connect very well as seen in figs. EHi discussed below. The results also match a 3rd 
code based on finite differences and a Newton-Krylov iteration to steady stat^^ (not shown here). Mirror 
symmetry (11) is imposed and eliminates the mean flow (1^. Shift-reflect symmetry (12) is not imposed 


explicitly but is satisfied by the steady solutions presented here. 

We show two branches of nonlinear steady solutions that bifurcate from the conduction state at Ra « 1708, 
a « 1.558. The primary branch has fixed horizontal wavenumber a (rounded here to a = 7r/2 ^ L/H = 2). 
The optimum branch adjusts a to maximize heat flux Nu. The (Ra, Nu) curves for both solutions are shown 
in figures and The steady primary solution is stable up to Ra « 53 000 where it spawns a time-periodic 
solution. The maximum and minimum Nu achieved by that periodic solution are shown in fig. that time 
periodic solution was computed with the time-accurate code 1. The unstable steady state was continued 
with code 2 up to Ra = 5.10® (symbols o in fig. [^. Mirror symmetry 0 suffices to stabilize the optimum 
solution in its fundamental periodic domain. That optimum solution was calculated up to Ra = 1.4 10® with 
code 1 and Ra = 10® with code 2 (symbols * in fig. 0. The optimum solution at Ra = 10® is well-resolved 
with 200 Chebyshev polynomials in y and 200 Fourier modes in x, after dealiasing. 

The heat flux for the unstable primary solution scales as Nu — 1 « 0.143 according to a least square 

fit in 5.10® < Ra < 5.10®. That solution develops a delicate spiral structure in the temperature field but 
not in the velocity (fig. . The winding of these temperature spirals continuously increases with Rayleigh 
number. Convergence of our algorithm to that unstable steady solution becomes quite slow for increasing 
Ra, apparently because of the center region of these spiral structures where v,T kQ. 

The optimum solution increases wavenumber a with Ra as a « 0.133 (fig.[^ to achieve an optimum 

heat flux Nu — 1 « 0.115 (fig. [^, from a least square fit in 10^ < Ra < 10®. This optimum heat 

flux scaling is quite similar to the heat flux observed in 3D turbulent convection experiments. Indeed, the 
optimum branch in fig. [ 2 ] is only slightly above the cryogenic helium gas data Nu « 0.088 RaP'^^ and 
the pressurized SFq gas data filfl^iVu « 0.105 (both in a cylinder of aspect ratio 1/2 and both for 

Ra < 10^^). The aspect ratio 4 datcP^ (Table 1) is well-fitted by Nu — 1 « 0.102 in 10® < Ra < 10^® 

and lies just below the optimum transport data with Nu — 1 « 0.115 This close agreement is 

remarkable given the differences in dimension, 2D optimum vs. 3D data, and Prandtl number, Pr = 7 for 
the 2D optimum vs. Pr « 0.7 in the helium gas experiments. 
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FIG. 2. Nu — 1 vs. Ra for primary branch {a = 7r/2, red o’s) with Nu — 1 ~ 0.143 (red dash) and optimum 

branch (red *’s) with Nu — 1 ~ 0.115(least square ht in 10^ < Ra < 10®, red solid). Lower (green) dash is 
the 3D turbulent data gjm jYu ~ 0.088 , (green) dash-dot is the 3D turbulent data filP^Afit ~ 0.105 

both for domain aspect ratio 1/2. The blue D’s for Ra > 10® is the aspect ratio 4 datj^ (Table 1, Nucorr)- 
Line ‘fs’ is the best free-slip upper bouncpSl Nu — 1 < 0.106 Line ‘ns’ is the best no-slip upper bouncP^ 

fVu - 1 < 0.02634 



FIG. 3. Horizontal wavenumber a = a°^^{Ra) that maximizes heat flux Nu, « 0.133 i?a®'®^^(least square fit in 
10'^ <Ra< 10®). 

The optimizing wavenumber a = a°P^{Ra) (fig.|^ shows an undulation between 1708 < Ra < 10® that is 
linked to a temperature spiral structure in the optimum transport solution as well (fig. [^. A temperature 
updraft of hot fluid (and downdraft of cold fluid) develops between Ra = 1708 and Ra « 10^, but a 
spiral structure with slight downdraft of warm fluid (and updraft of cool) develops in 10^ < Ra < 10®, 
corresponding to the bump in the curve a = a°P*{Ra) in fig. The warm spiral begins winding back up 
at Ra « 10®. The spiral does not appear to ever wind back down for higher Ra, developing instead an 
increasingly jagged structure (fig. right), and a°P* approaches the power law scaling 0 °^* ~ 0.133 
The structure of optimum Nu solutions depends on Prandtl number Pr. This is investigated in forthcoming 
worlP^ which also shows that the locally optimum steady solution discussed here is in fact the global optimum 
over a. 

The close agreement between the optimum transport 2D steady state solutions and 3D turbulent data in 
fig-i is intriguing. Although this agreement could be fortuitous, it strongly suggests that a single unstable 
steady solution may capture key statistical features of fully developed turbulent flows, such as the net heat 
flux and t he mean temperature profile as well as the strength and scale of elemental plumes, as in shear 
flow^HHni A search for such maximum momentum transport solutions in shear flows was initiated over 10 
years ago but not completed because of the higher computational complexity of those 3D solutions (J. Wang 
and F. Waleffe, 2004, unpublished). 

Our results are connected with Malkus’ theory of turbulent convectioii^ and subsequent work on upper 
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FIG. 4. Primary solution. Temperature T (top) and vertical velocity v (bottom) at Ra = 5. 10® for L/H = 2, 
Nu = 11.93. Equispaced contours at 10% of max-min, actual aspect ratio with —2 < x < 2 horizontal and 
< 2/ < 1 vertical. 



FIG. 5. Optimum solution. Velocity v (left) and temperature T (center) at Ra = 5 10®, L/H = 0.803, Nu = 14.72. 
Temperature T (right) at Ra = 10®, L/H = 0.262, Nu = 72.3. Equispaced contours at 10% of max-min, actual 
aspect ratio. 


Two key ingredients of Malkns’ theory are maximum heat transport and marginal stability 
of the mean temperature profile and the smallest scales of motion. Both ingredients are included in our 
calculations that maximize heat transport over horizontal wavenumber a and track steady state solutions 
that bifurcate from the marginal stability critical point at Ra ~ 1708, a ~ 1.558. We conjecture that our 2D 
results are in fact the 3D optimum transport solutions (for infinite or periodic horizontal directions). The 
upper bound results also assume 2D optimizers. Whether the ultimate scaling of these optimum Boussinesq 
solutions is Nu ~ as Ra —>■ oo remains to be seen but is possible, as is an abrupt transition to a 

smaller scale optimum solution.l^ 

Why the optimum transport solution should capture gross 3D turbulence characteristics might be under¬ 
stood as a ‘winner-take-air effect, where the optimum solution consumes all available potential energy so 
no other flow can be sustained. The optimum solution appears to be stable when mirror symmetry 0 
is imposed and length scales are restricted to be less or equal to the optimum wavelength. The optimum 
solution is unstable to larger scale perturbations, in particular to subharmonics where plumes merge and 
form bigger plumes in a cyclic or quasi-cyclic fashion. It is also unstable to a mean shear flow Q when 
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mirror symmetry is allowed to be broken. These instabilities would be the source of the Turbulence’ but 
the underlying unstable coherent solutions control the heat transport. 
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